/*
*		
*		A compilacao e feita por: gcc -o borwein_pthread borwein_pthread.c -lgmp -lm -lpthread
*
*/

#include <stdio.h>
#include <time.h>
#include <limits.h>
#include <gmp.h> // gcc -lgmp
#include <pthread.h>

#define N 25
#define BASE 10
#define MANTISSA 1024*1024*32

#define N_THREADS (N*4)


// define as funcoes
void *a_zero(); 
void *y_zero(); 
void *eq_exp(void *n);
void *eq_ay4();
void *eq_yep();
void *y_kn(); 

mpf_t expt;
char isExpInit; // Serve para verificar se exp ja foi inicializado

mpf_t ay4;
char isAy4Init; // Serve para verificar se ay4 ja foi inicializado

mpf_t yep;
char isYepInit; // Serve para verificar se yep ja foi inicializado

// define as variaveis
mpf_t a;
mpf_t y;


int main(int argc, char* argv[])
{
	unsigned long int k = 0, i = 0, t = 0;
	time_t t_zero, t_final;
	double delta_t;
    
	// inicializa a variavel
	mpf_t pi;

    pthread_t threads[N_THREADS];
    pthread_attr_t attr;

	// marca o tempo
	time(&t_zero);

	// Colocar Pthreads como Joinable    
    pthread_attr_init(&attr);
    pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);

	 // Float com precisao arbitraria
	mpf_set_default_prec(MANTISSA);

    isExpInit=0;
    isAy4Init=0;
    isYepInit=0;
	
    // mpf_t a && mpf_t y são inicializdos dentro das threads
    pthread_create(&threads[i], &attr, a_zero, NULL); i++;
    pthread_create(&threads[i], &attr, y_zero, NULL); i++;

	// sincroniza
    pthread_join(threads[t], NULL); t++;
    pthread_join(threads[t], NULL); t++;
    
	/*

	a = 6 - (4 * 2^1/2)
	y = (2^1/2) - 1

	yn = (1-(1-y^4)^1/4)/
	     (1+(1-y^4)^1/4)

	an = a(1+ yn)^4 - 2^(2k+3) * yn(1 + yn + yn^2)
	**	an converge para 1/pi

	*/
	
	for(k = 0; k < N; k++)
	{
		// cria a thread y_kn
        pthread_create(&threads[i], &attr, y_kn, NULL); i++;
		// cria a thread eq_exp
        pthread_create(&threads[i], &attr, eq_exp, (void*)k); i++;

        pthread_join(threads[t], NULL); t++; // espera y_kn
        
		// cria a thread eq_ay4
        pthread_create(&threads[i], &attr, eq_ay4, NULL); i++;
 
        pthread_join(threads[t], NULL); t++; // espera eq_exp
 
		// cria a thread eq_yep
        pthread_create(&threads[i], &attr, eq_yep, NULL); i++;

        pthread_join(threads[t], NULL); t++; // espera q_ay4
        pthread_join(threads[t], NULL); t++; // espera eq_yep

		// a = ay4 - yep
        mpf_sub(a, ay4, yep);    
    }

	// libera a memoria
	mpf_clear(y);
	mpf_clear(yep);
	mpf_clear(ay4);
	mpf_clear(expt);

    pthread_attr_destroy(&attr);
	
	// an converge para 1/pi
	mpf_ui_div(a, 1, a); 
	
	//mostra o valor de a
	mpf_out_str(NULL, BASE, 0, a);
	// libera a memoria
	mpf_clear(a);
	
	// mostra o tempo
	time(&t_final);
	delta_t = difftime(t_final, t_zero);
	printf("\nTempo = %.2lf segundos\n", delta_t);	
	
	// espera todas as threads acabarem de rodar antes da main terminar
    pthread_exit(NULL);

	return 0;
}


// calcula a_zero e recebe a variavel 'rop' por parametro
// a_zero = 6.0-(4.0*sqrt(2));
void a_zero(mpf_t rop)
{ 
	mpf_t aux;
	mpf_init(aux);
	
	// aux = sqrt(2)
	mpf_sqrt_ui(aux, 2); 
	// aux = aux * 4
	mpf_mul_ui(aux, aux, 4);
	// rop = 6 - aux
	mpf_ui_sub(rop, 6, aux);
	
	// libera a memoria
	mpf_clear(aux);
}


// calcula y_zero e recebe a variavel 'rop' por parametro
// y_zero = sqrt(2)-1.0;
void y_zero(mpf_t rop)
{ 
	mpf_t aux;
	mpf_init(aux);
	
	//aux = sqrt(2)
	mpf_sqrt_ui(aux, 2); 
	// rop = aux - 1
	mpf_sub_ui(rop, aux, 1);

	// libera a memoria
	mpf_clear(aux);
}


// calcula 2^(2k+3) e recebe a variavel 'n' por parametro
void *eq_exp(void *n)
{
    unsigned long int k;

    if(!isExpInit)
	{
        mpf_init(expt);
        isExpInit=1;
    }

    k = (unsigned long int) n;
     
	// k = (2k+3)
    k = (2*k)+3;
	// expt = 2
	mpf_set_ui(expt, 2);
	// expet = expt ^ k
	mpf_pow_ui(expt, expt, k);	

    pthread_exit(NULL);
}


// calcula y_kn e recebe as variaveis 'rop' e y por parametro
//	yn = (1-(1-y^4)^1/4)/
//	     (1+(1-y^4)^1/4)
void y_kn(mpf_t rop, mpf_t y)
{ 
	mpf_t y4, n, d;	
	mpf_init(y4);	
	mpf_init(n);
	mpf_init(d);
	
	
	// y4 = y^4
	mpf_pow_ui(y4, y, 4);
	// y4 = 1 - y4
	mpf_ui_sub(y4, 1, y4);
	// y4 = sqrt(y4)
	mpf_sqrt(y4, y4);
	// y4 = sqrt(y4)
	mpf_sqrt(y4, y4);
	
	// n = 1 - y4
	mpf_ui_sub(n, 1, y4);

	// d = 1 + y4
	mpf_add_ui(d, y4, 1);

	// rop = n / d
	mpf_div(rop, n, d);

	// libera a memoria
	mpf_clear(y4);
	mpf_clear(n);
	mpf_clear(d);
}


// calcula a(1+ yn)^4 
void *eq_ay4()
{
    if(!isAy4Init){
        mpf_init(ay4);
        isAy4Init=1;
    }
 
	// ay4 = y + 1
	mpf_add_ui(ay4, y, 1);
	// ay4 = ay4^4
	mpf_pow_ui(ay4, ay4, 4);
	// ay4 = a * ay4
	mpf_mul(ay4, a, ay4);

    pthread_exit(NULL);
}


// calcula 2^(2k+3) *  yn(1 + yn + yn^2)
void *eq_yep()
{
    mpf_t yk;

    if(!isYepInit)
	{
        mpf_init(yep);
        isYepInit=1;
    }
 

	mpf_init(yk);
	// yk = y^2
	mpf_pow_ui(yk, y, 2);
	// yk = yk + y
	mpf_add(yk, yk, y);
	// yk = yk + 1
	mpf_add_ui(yk, yk, 1);

    // yk = y(yk)
	mpf_mul(yk, yk, y);	

	// yep = expt * yk
	mpf_mul(yep, expt, yk);

	// libera a memoria
    mpf_clear(yk);

    pthread_exit(NULL);
}
